% Taylor diagram for summarizing model performance
%
% ---- References:
% Taylor, 2001 - Journal of Geophysical Research, 106(D7).
% Charria et al., 2008 - Biogeosciences, 5, 1437-1455.
% Jolliff et al., 2009 - Journal of Marine Systems, 76, 64-82.
%
% Program adapted from the IDL routine from K.E. Taylor
%
% ---- Call :
%    taydiag(tsig,rsig,tcorr,out,name_experiment,title)
%
% ---- Input:
%    Needed:
%       tsig     : Standard deviation for the experiments (n values)
%       rsig     : Standard deviation for the reference (1 or n values)
%       tcorr    : Correlation between the reference and the experiments
%                 (n values)
%
%   List of options:
%   - 'MarkerDisplayed' : 'marker' (default) / 'colorbar'
%       -> marker: Experiments are represented by individual symbols
%       -> colorbar: Experiments are represented by a color described in a
%       colorbar
%
%   - 'z' : When 'MarkerDisplayed'='colorbar', another dimension can be
%           described (e.g. depth or longitude). The vector must have n
%           values as tsig and tcorr.
%
%   - 'colorbartitle' : Title for the colorbar. By default, empty string.
%           This option is ignored if 'MarkerDisplayed' == 'marker'
%
%   - 'npan' : 1 or 2 panels (Quarter-circle) to display
%       (1 for positive correlations,
%				2 for positive and negative correlations)
%				Default value depends on tcorr values.
%
%   - 'rmslines' : 'off' (default) / 'on'
%       -> When 'on', plot the unbiased RMS difference lines
%
%   - 'out' : 'window' (default) /'file'
%       -> window: standard display in a window (default)
%       -> file: encapsulated postcript output named 'output.eps'
%
%   - 'filename' : EPS file name when printed output. Default file name is
%   'output.eps'
%
%   - 'name_experiment' : Name of the experiment used in the legend
%                 ([n x length of the longer name] array of characters)
%                       -- Or cell with experiment names
%                           >> name_experiment={'BATS','EUMELI'};
%
%   - 'title' -> Title to add (string)
%
%   - 'arccntl' : 0 / 1 (default)
%       -> If arccntl == 1, radial lines are plottes
%       -> If arccntl == 0, only ticks are visible
%
% ---- Output:
%       postcript or figure depending on the 'out' option ...
%
%
% ----------------------------------------------------------------------
% Created: G. Charria (02/2008)
% -- Updated: 03/2010 --
% v0 : Code d'origine
% v1.0 : - polar coordinates (03/2010)
%        - option 'filename'
%        - 180° or 90° plot (for negative correlations)
%        - option 'arccntl'
%        - option 'MarkerDisplayed' + 'z'
%        - improvement of legend position
%        - option 'rmslines'
% ----------------------------------------------------------------------
function taydiag(tsig,rsig,tcorr,varargin)

%close all

% ---------------------------
% ---- Check input
% ---------------------------
if nargin < 3, help(mfilename), return, end
narg=nargin - 3;
if mod(narg,2) ~=0
    error('Wrong number of input arguments (see help taydiag) !')
end
% -- Scan options
% - Default parameters
out='window';
maintit=[];
name_experiment=[];
filename='output.eps';
index=find(tcorr<0, 1); % -- Control the sign of correlations
if ~isempty(index)
    Npan=2;
else
    Npan=1;
end
arccntl=1;
markerdisplayed='marker';
colorbartitle='';
z=[];
rmslines='off';

for iopt = 1 : 2 : narg
    optname  = varargin{iopt};
    optvalue = varargin{iopt+1};
    switch lower(optname)
        case 'out'
            out=optvalue;
        case 'title'
            maintit=optvalue;
        case 'filename'
            filename=optvalue;
        case 'npan'
            Npan=optvalue;
        case 'name_experiment'
            name_experiment=optvalue;
        case 'arccntl'
            arccntl=optvalue;
        case 'markerdisplayed'
            markerdisplayed=optvalue;
        case 'z'
            z=optvalue;
        case 'colorbartitle'
            colorbartitle=optvalue;
        case 'rmslines'
            rmslines=optvalue;
    end
end

% -- Control on the reference rsig
index=find(rsig==0);
if ~isempty(index)
    rsig(index)=NaN;
    disp(['Warning: ',num2str(length(index)),' "rsig" values are equal to zero and are not displayed !'])
end

% ---------------------------
% ---- Parameters definition
% ---------------------------
% sigmax : radial limit of plot
if (max(tsig./rsig) < 1)
    sigmax = 1;
else
    sigmax = max((tsig./rsig))+0.1;
end
% minsigmax : left limit of plot
% mincorr : minimum correlation coefficient
if Npan==1
    minsigmax=0;
    mincorr=0;
else
    minsigmax=-sigmax;
    mincorr=-1;
end


% Y-axis title
ytitle = 'Standard Deviation (Normalized)';
% X-axis title
if Npan==1
    xtitle = '';
elseif Npan==2
    xtitle = 'Standard Deviation (Normalized)';
end
% Fontsize for the plot
axcsiz = 12;
% Tick location
yticloc=0.:.1:sigmax;
ylabels = yticloc;
if Npan==1
    xticloc = yticloc;
    xlabels = ylabels;
elseif Npan==2
    xticloc = [-yticloc(end:-2:2),yticloc(1:2:end)];
    xlabels = xticloc;
end

dummy=0.;

% -------------------------------------
% Ticks
% -------------------------------------
switch Npan
    case 1
        % Minor Ticks
        rminticval = 0.05:0.1:0.95;
        % Subminor Ticks
        rsubticval = 0.91:0.01:0.99;
        % Major Ticks
        rlabticval = [0.0, 0.1, 0.2, 0.3, 0.4,...
            0.5, 0.6, 0.7, ...
            0.8, 0.9, 0.95, 0.99, 1.0];
        % Major Tick Labels
        rticlabels = {'0.0', '0.1','0.2',...
            '0.3', '0.4', '0.5',...
            '0.6', '0.7', '0.8',...
            '0.9', ' 0.95', ' 0.99',...
            '1.0'};
    case 2
        % Minor Ticks
        rminticval = -0.95:0.1:0.95;
        % Subminor Ticks
        rsubticval = [-0.99:0.01:-0.91 0.91:0.01:0.99];
        % Major Ticks
        rlabticval = [-1.0, -0.99, -0.95, -0.9, -0.8,...
            -0.7, -0.6, -0.5,...
            -0.4, -0.3, -0.2, -0.1,...
            0.0, 0.1, 0.2, 0.3, 0.4,...
            0.5, 0.6, 0.7, ...
            0.8, 0.9, 0.95, 0.99, 1.0];
        % Major Tick Labels
        rticlabels = {'-1.0','-0.99','-0.95','-0.9',...
            '-0.8','','-0.6',...
            '','-0.4','','-0.2','',...
            '0.0','','0.2',...
            '','0.4','',...
            '0.6','','0.8',...
            '0.9','0.95','0.99',...
            '1.0'};
end

% 'Correlation' label
corrt=['C';'o';'r';'r';'e';'l';'a';'t';'i';'o';'n'];
% Length of minor ticks
minortic= 0.015;
% Length of subminor ticks
subtic = 0.010;

% -------------------------------------
% Definition of markers
% -------------------------------------
switch markerdisplayed
    case 'marker'
        kind=['+';'o';'x';'s';'d';'^';'v';'p';'h';'*'];
        colorm=['b';'r';'g';'c';'m';'y';'k'];
        n=1;
        marker(1:(size(colorm,1)*size(kind,1)),1:2)=' ';
        for ic=1:size(colorm,1)
            for ik=1:size(kind,1)
                marker(n,:)=[kind(ik,:) colorm(ic,:)];
                n=n+1;
            end
        end
        if (length(tsig) > 70)
            disp('You must introduce new markers to plot more than 70 cases.')
            disp('The ''marker'' character array need to be extended inside the code.')
            return
        end
    case 'colorbar'
end

% --------------------------------------
% Test to draw main rays or not
if (arccntl > 0)
    rticlen=1.; % Draw the complete ray
else
    rticlen=0.02; % Draw just a Tick
end



figure;
switch out
    case 'file'
        set(gcf,'Visible','off')
end

% -----------------------------
% ---- Plot main axes (X and Y)
% -----------------------------
plot(dummy,dummy)
box('off')

if Npan==1
    pbaspect([1 1 1])
else
    pbaspect([2 1 1])
    % __ Remove default axes
    set(gca,'YColor','w',...
        'YTick',[],...
        'YTickLabel',[])
end


set(gcf,'Color','w')
set(gca,'XLim',[minsigmax sigmax],...
    'YLim',[0., sigmax],...
    'YMinorTick','on',...
    'XMinorTick','on',...
    'TickLen',[0.01 0.01],...
    'XTick',xticloc,...
    'XTickLabel',xlabels,...
    'YTick',yticloc,...
    'YTickLabel',ylabels,...
    'FontSize',axcsiz)
xlabel(xtitle)
ylabel(ytitle)

hold on;

% -----------------------------
% ---- Plot radial lines
% -----------------------------
% Plot the extreme line
[xa,ya]=pol2cart(acos(mincorr:0.001:1),sigmax);
plot(xa,ya,'k')

% Plot the '1' line
[xa,ya]=pol2cart(acos(mincorr:0.001:1),1.);
plot(xa,ya,':k','LineWidth',1.5)

% Plot few other lines ...
for il=0.2:0.2:1
    [xa,ya]=pol2cart(acos(mincorr:0.001:1),il);
    plot(xa,ya,':k')
end
for il=1:(sigmax-1)/3:sigmax
    [xa,ya]=pol2cart(acos(mincorr:0.001:1),il);
    plot(xa,ya,':k')
end

% Plot the Y centered axis
if Npan==2
    plot([0 0],[0 sigmax],'k')
end

% -----------------------------
% ---- Plot RMS lines
% -----------------------------
if strcmp(rmslines,'on')
    rms_value=0:.2:abs(1-minsigmax);
    for irms=1:length(rms_value)             
        [xrms,yrms]=pol2cart(acos(-1:0.01:1),rms_value(irms));   
        xrms=xrms+(rsig(1)/rsig(1));
        [th,rho]=cart2pol(xrms,yrms);
        rho(rho>sigmax)=NaN;
        [xrms,yrms]=pol2cart(th,rho);
        plot(xrms,yrms,'LineStyle','-.','Color',[0 0.5 0])
        % Attention: normalement 1 en normalise mais rsig pour preparer le
        % cas non normalise.
        
        % ---- Ticks RMS
        if irms > 1
            text(xrms(3),yrms(3),num2str(rms_value(irms)),...
                'BackgroundColor','w','Color',[0 0.5 0],...
                'HorizontalAlignment','center',...
                'Rotation',75);
        end
    end
end


% ------------------------------------
% ---- Draw tick marks (on the circle)
% ------------------------------------
% *** Main Ticks *** + Correlation labels + 'Correlation' title
for i=1:length(rlabticval)
    
    r = rlabticval(i);
    [ticx(1),ticy(1)]=pol2cart(acos(r),sigmax);
    [ticx(2),ticy(2)]=pol2cart(acos(r),(sigmax*(1.-rticlen)));
    
    if i==1 || i==length(rlabticval)
        plot(ticx,ticy,'k');
    else
        h=plot(ticx,ticy,'b');
        set(h,'LineStyle','-.')
    end
    
    % Tick Labels
    if Npan==1
        [xt,yt]=pol2cart(acos(r),(sigmax+0.04));
    elseif Npan==2
        [xt,yt]=pol2cart(acos(r),(sigmax+0.1));
    end
    tt=text(xt,yt,rticlabels(i),...
        'Fontsize',axcsiz,...
        'Color','b',...
        'HorizontalAlignment','center');
    if rlabticval(i)==0
        b=get(tt,'Position');
        sfly=b(2)+0.03;
    end
    if rlabticval(i)==1
        a=get(tt,'Parent');
        b=get(a,'Position');
        sfl=b(1)+b(3)+0.03;
    end
end

% *** 'Correlation' title ***
for i=1:size(corrt,1)
    r=0.7+(i-1)*((0.9-0.7)/12);
    if Npan==1
        [xt,yt]=pol2cart(acos(r),(sigmax+0.1));
    elseif Npan==2
        [xt,yt]=pol2cart(acos(r),(sigmax+0.2));
    end
    text(xt,yt,corrt(i,:),...
        'Fontsize',axcsiz,...
        'Rotation',-((pi/2)-acos(r))*180/pi,...
        'Color','b')
end

% *** Minor Ticks ***
for i=1:length(rminticval)
    
    r = rminticval(i);
    [ticx(1),ticy(1)]=pol2cart(acos(r),sigmax);
    [ticx(2),ticy(2)]=pol2cart(acos(r),(sigmax*(1.-minortic)));
    plot(ticx,ticy,'k')
end

% *** Sub minor Ticks ***
for i=1:length(rsubticval)
    r = rsubticval(i);
    [ticx(1),ticy(1)]=pol2cart(acos(r),sigmax);
    [ticx(2),ticy(2)]=pol2cart(acos(r),(sigmax*(1.-subtic)));
    plot(ticx,ticy,'k')
end

% ------------------------------------
% ---- Draw Results
% ------------------------------------
% Loop over the number of experiments
hp(1:length(tsig))=0.;
for i1=1:length(tsig)
    if length(rsig) > 1
        
        % Ancienne formulation ... erreur negligeable (O(1e-16)) par rapport a polaire
        %xxcontrol = (tsig(i1)/rsig(i1))*tcorr(i1);
        %yycontrol = (tsig(i1)/rsig(i1))*sqrt(1.-tcorr(i1)*tcorr(i1));
        
        % cf. Figure 1 in Taylor paper
        [xx,yy]=pol2cart(acos(tcorr(i1)),(tsig(i1)/rsig(i1)));
        
    else
        
        % cf. Figure 1 in Taylor paper
        [xx,yy]=pol2cart(acos(tcorr(i1)),(tsig(i1)/rsig(1)));
        
    end
    
    switch markerdisplayed
        case 'marker'
            hp(i1)=plot(xx,yy,marker(i1,:),'MarkerSize',8,...
                'MarkerFaceColor',marker(i1,2),...
                'Linewidth',2.5);
            if (strcmp(marker(i1,1),'o')||...
                    strcmp(marker(i1,1),'s')||...
                    strcmp(marker(i1,1),'d')||...
                    strcmp(marker(i1,1),'^')||...
                    strcmp(marker(i1,1),'v'))
                set(hp(i1),'MarkerEdgeColor','k')
            end
            %set(hp(i1),'DisplayName',num2str(mC(1,k))
            
        case 'colorbar'
            if ~isempty(z)
                hp=scatter(xx,yy,50,z(i1),'filled');
            else
                hp=scatter(xx,yy,50,i1,'filled');
            end
            
            
    end % end switch markerdisplayed
    
    
end


% ------------------------------------
% ---- Draw Legend / Colorbar
% ------------------------------------
% -- Colorbar
if strcmp(markerdisplayed,'colorbar')
    hc=colorbar('EastOutside');
    posc=get(hc,'Position');
    outposc=get(hc,'OuterPosition');
    tightposc=get(hc,'TightInset');
    annot=get(hc,'Title');
    set(annot,'String',colorbartitle)
    set(hc,'Title',annot)
    set(hc,'Position',[posc(1)-outposc(1)+tightposc(1) posc(2)-outposc(2)+tightposc(2) posc(3)/2 posc(4)])
end
% -- Legend
overx=0;
if ~isempty(name_experiment) && strcmp(markerdisplayed,'marker')
    
    ll=legend(hp,name_experiment,'location','EastOutside');
    pleg=get(ll,'Position');
    overx=sfl+pleg(3)-1;
    if overx <= 0
        overx=0;
    end
    set(ll,'Position',[sfl-overx*sfl pleg(2) pleg(3) pleg(4)])
end
% ------------------------------------
% ---- Title
% ------------------------------------
if ~isempty(maintit)
    htit=title(maintit,'Fontsize',axcsiz);
    ptit=get(htit,'Position');
    set(htit,'Position',[ptit(1) sfly])
end

% ------------------------------------
% ---- Resize step if needed
% ------------------------------------
fpos=get(gca,'Position');
fpos(3)=fpos(3)-overx*fpos(3);
fpos(4)=fpos(4)-overx*fpos(4);
set(gca,'Position',fpos)
%pleg(3)/

% ------------------------------------
% ---- Print eps file
% ------------------------------------
switch out
    case 'file'
        set(gcf,'PaperOrientation','portrait')
        set(gcf,'PaperUnits','centimeters')
        set(gcf,'PaperType','A4')
        % Modification to display the whole legend
        if strcmp(markerdisplayed,'marker')
            set(gcf,'PaperPosition',[1. 1. 38.+(10*(max((length(tcorr)/10),3) - 3)) 38.+(10*(max((length(tcorr)/10),3) - 3))])
        else
            set(gcf,'PaperPosition',[1. 1. 38. 38.])
        end
        print('-f','-depsc2','-r300',filename)
end
